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Abstract 

We introduce a time-integrator to sample with high order of accuracy the invariant 
distribution for a class of semilinear SPDEs driven by an additive space-time noise. 
Combined with a postprocessor, the new method is a modification with negligible over¬ 
head of the standard linearized implicit Euler-Maruyama method. We first provide an 
analysis of the integrator when applied for SDEs (finite dimension), where we prove that 
the method has order 2 for the approximation of the invariant distribution, instead of 1. 
We then perform a stability analysis of the integrator in the semilinear SPDE context, 
and we prove in a linear case that a higher order of convergence is achieved. Numerical 
experiments, including the semilinear heat equation driven by space-time white noise, 
confirm the theoretical findings and illustrate the efficiency of the approach. 

Keywords: stochastic partial differential equations, postprocessor, invariant measure, 
ergodicity, space-time white noise. 
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I Introduction 

We introduce an efficient integrator for the sampling of the invariant probability distribution 
of a class of semilinear parabolic SPDEs with additive noise written as an abstract stochastic 
evolution equation (in the sense of m) 

du{t) = {Au{t) + F{u{t))) dt + dW^{t), u{0) = uq. (1) 

Its solution u{t) takes values in a separable inhnite dimensional Hilbert space T-L, with 
initial condition uq (assumed deterministic for simplicity). We assume that —A is a positive 
unbounded self-adjoint linear operator with an associated sequence of positive eigenvalues 
0 < Ai < A 2 < ... and an associated complete orthonormal family of eigenvectors ei, 62 ,... 
The coefficient F : T-L ^ T-L is a Lipschitz continuous nonlinearity, and it is assumed to 
derive from a continuously differentiable potential function V : Ti ^ i.e. F = —DV. 
Finally, we assume that is a Q-Wiener process on Ti dehned on a probability 

space (D,J-', P) fulfilling the usual conditions; the covariance operator Q ■. Ti ^ Ti is a 
bounded, non-negative self-adjoint linear operator, such that Qei = qiCi for some bounded 
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sequence of real numbers {qi)i£n *, where we use the notation N* = {1, 2, 3,..We assume 
the following trace condition: 

s = sup |s G (0,1) : Trace< +oo| > 0. (2) 

Then, there exists a unique mild solution of ([T]) on R+ (see HD Chapter 7]), i.e. an'H-valued 
continuous stochastic process which satisfies 

u{t) = e^^uo+ [ e^^-^^^F{u{s))ds+ [ e^^-^^^dW^{s). (3) 

Jo Jo 

This abstract setting includes the stochastic semilinear heat equation where A = A = 
in Ti = L'^i'D) for some open smooth bounded domain 7) C 


du, . . , , , ,, dW'^ . . 

— (a:,t) = Au{x,t) + f{u{x,t)) + —^{x,t), 


( 4 ) 


with homogeneous Dirichlet boundary conditions on (9P. In ([3|), the coefficient / ; M —>■ M 
is a smooth, Lipschitz function; the associated Nemytskii coefficient F: u^T-Lt-^fou&'H 
satisfies F{u) = —DV{u) for all u with V{u) = — Jq {F{9u), vJ'^dO. 

If d = 1, one can consider space-time white noise in ([3]), (he. with the identity covariance 
operator Q = F, this choice yields s = 1/2 in (l2|)). When d > 1, a nontrivial covariance 
operator Q 7 ^ / is required, yielding noise which is white in time and colored in space. Notice 
also that (j2]) is automatically satisfied with s = 1 for a trace-class noise with Trace(Q) < 00. 

To study the long-time behavior of the process u, we make the additional assumption 
that T in ([1]) is Lipschitz continuous with constant L > 0 such that L < Ai = miupgf^* Ap. 
Then (see e.g. Section 8.6 in [TU]) Eq. (pQ) admits a unique invariant distribution //oo- This 
means that for all (smooth and Lipschitz) test functions 1 /) : "H —)■ M, and for all initial 
conditions uq , 

lima.s. (j){u{t))dt = [ (j){y)di^^{y), 

T^oo T Jq 

where the notation lima. s. means that the limit holds with probability 1. Moreover (see 
e.g. Section 6.3 in m ) u{t) converges in law to /ioo exponentially fast in the following 
sense: for all t > 0, and any test function cj), 


K(j){u{t)) 



(l){y)dy 

00 (d) 


<C'((/,uo)e 


( 5 ) 


where C{4>,uo) is independent of t. It is a standard approach to take advantage of the 
estimate ([S]) to compute ergodic integrals of the form (/)(y)d/ioo(y)- To do so, in practice 
one needs to rely on a discretization of the evolution equation. We thus now present two 
implicit-explicit time-discretization schemes. 


Linearized implicit Euler method We first consider the simplest numerical scheme, 
which is referred to as the linearized implicit Euler method in this article. Given a constant 
timestep size h > 0, it is defined by Vn+i = Vn + hAvn+i + hF{vn) + Vh^n , equivalently 

Vn+l = Jl (vn + hF{Vn) + , (6) 
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where vq = u(0) = uq, Ji = (/ — hA) ^ and in = h + l)h) — W'^{nK)). For 

a fixed final time T > 0 , it is known that the scheme ([ 6 ]) applied to © has weak order of 
accuracy q for all q <s, i.e. it satisfies for all h small enough, 

|]E(,^(n„)) - E(<^(n(t„)))| < T)h\ (7) 

for all tn = nh < T, where C{uq, 4>, T) is independent of n, h. 

The weak order of accuracy d?]) has been analyzed in m for ([ 6 ]) applied to the heat 
equation in the case F = 0 and in |12] in the semilinear case; see also [39] for the case of 
colored noise in dimension d > 1 and |23] where different techniques are used. Concerning 
the approximation of the invariant distribution of the heat equation, it is proved in | 1 ] that 
for d = 1 and Q = I, one has for all time tn, similarly to dSj), the exponential convergence 
property 

E(</>(^n)) - / 4>iy)dtLoo{y) 

Jh 

for all r < s with s = 1 / 2 , where K{uo,(j)) are independent of n and h, and A is 

independent of (f), h, n. The proof is based on the analysis of the weak approximation using 
the tools of |12] (expansion of the error using the backward Kolmogorov equation, Malliavin 
calculus), with a proof that constants C{uo, (p,T) in ([7|) can be chosen independent of the 
final time T, in the spirit of |35| . 

New method The contribution of this paper is the introduction of a modification, de¬ 
noted Un, of the standard Euler scheme (l 6 |), together with a postprocessor, denoted Un, 
which permits to achieve the estimate dS]) with the higher order r < s -|- 1 instead of r < s, 
when in ([ 8 |) is replaced with Precisely, the postprocessed method is given by two 
sequences and in % defined by 

Un+l = Jl (un + hF{un + ^VhJ2in) + ^ ^ (9) 

Un = Un + ( 10 ) 

where we introduce additional operators J 2 , J 3 as follows: 

J2 = (/ - ^-^hA)-\ JsQJi = {I- ^A)-^Q. 

Note that J 3 is not determined uniquely by the equation above: one may define J 3 = 
(/ — or use a Cholesky decomposition. We emphasize the importance of the 

postprocessing operation (| 10 p : indeed, the order of convergence to the invariant distribution 
yLoo is not improved if in ([ 8 |) is replaced with Un, instead of Un- Notice that the new 
scheme ([^- (fTOjl has a negligible overcost compared to (| 6 |). First, the postprocessor (fTOjl 
needs not to be computed at each time step, but it can be computed only once at the end of 
each numerical trajectory. Second, computing J 2 in at each timestep is the only additional 
calculation. Third, we prove that the constant A > 0 in dHj) can be chosen with the same 
size for both methods: thus in terms of cost (number of timesteps required to achieve a 
given accuracy on the left-hand side of dSj)), our new integrator performs better than the 
standard scheme ([ 6 |). 


<K{uo,(P)e-^^- +Ci<PW ( 8 ) 
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The fact that the long-time weak order of accuracy r in ([5D can be made strictly larger 
than the short time accuracy g in ([7]) is not surprising: this is known for SDEs {i.e. in finite 
dimension in the terminology of this paper), see [2611241 fTl I38| in the context of Brownian 
dynamics and mm in the context of Langevin dynamics, where integrators with low weak 
order, typically g = 1 in ([7j), are shown to achieve a high order r > q over long times 
([ 8 |) for sampling the invariant measure. Inspired by these recent advances, the popular 
technique of processing for deterministic differential equations [ 8 ] was recently extended to 
the stochastic context in [38], and serves as a crucial ingredient to derive the new method 
proposed in this paper. In our context, this idea is to enhance the accuracy of the modified 
numerical method Un in ([9|) by applying a suitable change of variables Un defined in 

m- 

Alternatively, note that high order integrators in the strong sense (approximation of the 
trajectories instead of the distribution) for parabolic problems of the form ([1]) are proposed 
in |23[2niEI]; however these schemes belong to the class of exponential integrators, while 
the method proposed in this paper avoids the computation of matrix exponentials. We also 
mention another natural integrator for m, which is the (stochastic) trapezoidal method 
(also known as the Crank-Nicolson method), 

Un+l — Ufi —A.(^Un T -\- hF{un) T (f f) 

However, in contrast to ([U]), the scheme m is not L-stable, a desirable property for severely 
stiff problems, already in the deterministic literature [15] . This makes exponential conver¬ 
gence estimates of the form Q not true in general for the scheme (Hip . 

Finally, since the space T-L is infinite dimensional, a space discretization scheme is re¬ 
quired in practice, e.g. finite differences or finite elements - see Section 01 In this article, 
we only focus on the time-discretization issue. We mention m where a postprocessing 
technique is applied to improve the spatial discretization in the strong sense. 

Outline and main results This paper is organized as follows. In section [2l we explain 
the derivation of a generalized nonlinear version of the integrator (l^- dlOl) in the context of 
finite-dimensional SDEs, where many standard analysis tools are available compared to the 
SPDE context. We prove that the new integrator has order 2 for the approximation of the 
invariant distribution of nonlinear ergodic SDEs, instead of order 1 for the standard Euler 
scheme. Section 0] is devoted to the analysis of the integrator ([ ^ - (flUl) for the SPDEs ([I]): 
we detail the abstract Hilbert space setting ISection 13.111 . we show the stability and er- 
godicity properties of the integrator (Section 13. 2D . we prove in a simplified linear case the 
improved order of accuracy s-|- 1 of the new method for SPDEs, in contast to the order s for 
the standard Euler scheme (Section 13.3p . and we show that the proposed scheme exhibits 
the correct spatial regularity of the invariant distribution in terms of Sobolev-like spaces 
(Section 13.4p . Finally, Section 0] is dedicated to numerical experiments which confirm the 
theoretical findings and illustrate the efficiency of the new method. 
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2 New high order integrator: derivation and analysis in finite 
dimension 


In the context of SDEs in finite dimension N S N*, we consider the more general case of a 
nonlinear system of the form 


dX{t) = (/i(X(f)) + f2iXit)))dt + adW^it), X(0) = Xo, (12) 

with solution X[t) in . The nonlinearities /i ,/2 : —>■ are smooth and Lipschitz 

functions such that fi{x) + f 2 {x) = fo{x) = — VVb(x) for some potential function Vq : 

M, where /i is a term to be treated implicitly and /2 is a term to be treated 
explicitly. The initial condition Xq is assumed deterministic for simplicity. We also dehne 
the Q-Wiener process W^{t) = (t) where Q now denotes a NxN symmetric positive 

definite matrix and W{t) is a standard V-dimensional Wiener process, and u > 0 is a hxed 
constant. 

We introduce the following new implicit-explicit scheme for sampling with high order 
two the invariant measure of (m, 


Xn+l — Xn + hfi I Xn+l + 


-2 +V5 


Jr 


,2<TV^Cn ^ + hf2 (^Xn + ^^,20-V^Cn ^ 


, ^1-V2 + V5 ,_1 , 1 + V2-A . /T.O 

+ (- 2 - 2 - )Jn,2<XVh^^, 

— Xn + 1 

where J„,i, Jn,2, Jn,3 are given by 


(13) 


Jn,l = {I - hf[{Xn))-\ Jn,2 = {I- 






hf[{Xn))-\ Jn,3QJi3 = (^ " Ij/K^n))-^ 


We emphasize that the matrix inverses Jn,i, Jn, 2 , Jn ,3 are used only in the notations to 
define the scheme but should not be computed in practice. Indeed, in practical implementa¬ 
tions, a LU decomposition should be used in place of computing matrix inverses. Moreover, 
in the semilinear case (HI, this decomposition needs to be done only once and may be used 
for all further iterations. 


Remark 2.1. After a spatial discretization with finite differences or finite elements of the 
SPDE (JT]), in general one arrives at a system of stiff SDEs in with large dimension N 
of the form (fT^ , 

dX{t) = AX{t)dt + f{X{t))dt + adW^it). (14) 

where assume that f{x) = — W(a:) for some V : M. It corresponds to the special 

case /i(x) = Ax where —A now denotes a N x N symmetric positive definite matrix. In 
this case, the above scheme (USD simplifies to the integrator Q. Moreover, in the special 
case fi = 0, the method (HI reduces to 


Xn-{-l — Xji hf j Xrfi -\- j “1“ 


Xn= Xn + ^(rVh^n, (15) 


a method which was first proposed in and analyzed in IWf . although it was constructed 
in another manner using a non-Markovian formulation. 
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For simplicity, we assume for the remaining of this section that Q = I is the identity 
matrix, without loosing generality, applying the appropriate change of variable in (hxed) 
hnite dimension. Section O is devoted to a stability analysis in the case of an Ornstein- 
Uhlenbeck process: we prove L-stability and exactness results, which play an important 
role in the efficiency of the integrator for semilinear SPDEs. Section 12.21 contains details 
on the construction of the method in order to satisfy the above properties and to achieve 
order two of accuracy for sampling the invariant measure. 

2.1 Stability analysis 

For the study of the stability of stochastic integrators applied to stiff SDEs with additive 
noise of the form a widely used test problem is the following scalar SDE problem 

(Ornstein-Uhlenbeck process) 

dX = -XXdt + adW{t), (16) 

where A, a > 0 are fixed constants. Notice that this test equation provides only a useful 
insight but no rigorous general conclusion on the numerical long-time behaviour for nonlin¬ 
ear problems with additive noise, see [6] and references therein. A rigorous analysis of the 
proposed scheme in our semilinear SPDE context is presented in Section 13.21 Notice also 
that other test equations are used in the literature in the case of multiplicative noise, see 
[SlliaiMlElETlEaiS]. Consider a one step method of the form 

Xn+l = ■A{z)Xn + B{z)'/ha^n: z = “Ah, (17) 

where h is the stepsize, A{z),B{z) are analytic functions, ~ AA(0,1) are independent 
Gaussian random variables. The SDE (jl6l) is ergodic with the unique invariant measure 

a Gaussian with mean zero and variance i.e with density Pooix) = y^^^exp(— ^x^). 
Indeed, for any initial condition Xq = x, the exact solution is a Gaussian random variable, 
with lim^^oo E(-^(^)) = 0 and lim^^oo ]E(|Ai(t)p) = The second-order moment E(|X„p) 
remains bounded as n —)■ oo if |Al( 2 :)| < 1, which corresponds to the mean-square stability 
condition. The condition |AI(2:)| < 1 for all 2; with negative real part is called A-stability in 
the deterministic literature m- This is a desirable property of numerical integrators for 
stiff problems, because it permits to avoid a severe timestep size restriction. For |Al( 2 ;)| < 1 
(stability condition), we obtain 

lim E(|X„|2) = ^TZiz), where n{z) = (18) 

n^oo 2A 1 — A{z)^ 

We see that the method is exact (for the approximation of the invariant distribution) if 
and only if TZ{z) = 1 for all z] this is the case for instance for the trapezoidal method 
(fTTl) . which is such that A{z) = B(z) = However, in addition to A-stability, 

a desirable property of Runge-Kutta methods for very stiff problems is L-stability US], 
namely Al(oo) = 0; this is not satisfied by the trapezoidal method, for which Al(oo) = —1. 
The following proposition states that for Runge-Kutta type methods, where A{z),B{z) are 
rational functions, L-stability {i.e. Al(oo) = 0) is incompatible with the exactness for the 
invariant distribution, i.e. TZ{z) = 1. 

Proposition 2.2. Consider a method of the form (HZD where A{z),B{z) are rational func¬ 
tions. If the method samples exactly the invariant distribution of (fTHl) (i.e. IZ{z) = 1 then 
|.A(oo)| = 1. 
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Proof. If A{z) is a rational function with |^(oo)| A 1) there exist a constant C S M and an 
odd integer A: € Z such that TZ{z) ~ Cz^ for 2 ^ 00 , and thus 7l{z) ^1. □ 


However, as shown below, the above barrier for L-stable Runge-Kutta methods can be 
circumvented by applying an appropriate postprocessor to (I17p of the form 

Xn = C{z)Xn+ 'D{z)Vhaf,n, z = -Xh. (19) 

We will see that this feature serves as a crucial ingredient in Section [3] to achieve high order 
in the SPDE case. 

Proposition 2.3. The method (m applied to the SDE test problem (unp with fi{x) = 
—Ax and / 2 (x) = 0 is such that the scheme X^ e->■ Xn+i is L-stable and X^ is ex¬ 
act for sampling the invariant measure, i.e. for all test functions (p and all timesteps h, 

= lim„^+oo E((/>(W„)) = f^ 4 >{y)poo{y)dy. 

Proof. The method is L-stable because it has the same stability function A{z) = as 
the linearized implicit Euler method Q. A calculation using the above notations yields 

2 

lim E(X^) = 7r-r7?.(z), where TZ{z) = C{z)‘^TZ{z) — 2z'D{z)‘^, (20) 

n—)-cx> 2A 

with 2: = —Xh and TZ{z) given by (fTH]) . Noting that A{z) = B{z) = (1 — z)~^, C{z) = 
1, Viz) = ^(1 — zl 2 )~^A^ yields TZ{z) = 1, which proves that the method is exact. □ 


Remark 2.4. Relaxing the assumption that the matrices A and Q commute, the scheme 
(USD can he adapted to remain exact for sampling the invariant measure when applied to 
the linear problem dX = AXdt -\- aQ^/‘^dW without nonlinearity (fi{x) = Ax,f 2 {x) = 0). 
In fact, the definition of matrix J 3 should be modified, and obtained solving a system of 
Lyapunov eguations. 

2.2 Construction of the integrator (fT^ 

We explain in this section how we construct method (1131) with high order of accuracy for 
the invariant measure by using the idea of postprocessing from [38]. Consider a system of 
SDEs in of the form 

dX{t) = fo{X{t))dt + adW{t), A(0) = Xq, (21) 

where cr > 0 and {W(t))^^^_^_ is a standard d-dimensional Wiener process. 

Assumption 2.5. We assume the following. 

1. fo is of class C°°, with bounded derivatives of any order, and there exists a potential 

function Vq : —>• M such that fo = —VVo/ 

2. there exist C, /3 > 0 such that for all x € x'^/o(x) < —fix'^x C. 
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Then, system (l^TD is ergodic with a unique invariant measure /Xoo (see e.g. [TB]) given 


by 


/^oo(dy) = p{y)dy, with p{y) 


1 2 
exp(-^Vb( 2 /)). 


Z 


cr^ 


Moreover the solution X{t) of (I2ip satisfies an exponential ergodicity property analogous 
to (I5]). 

The following theorem is the main result in [35] where postprocessed integrators for 
SDEs are introduced. It permits to improve the accuracy of a method of weak order q 
to order r = q + 1 for the invariant measure. Notice that this theorem remains valid for 
general classes of SDEs with additive or multiplicative noise in multiple dimensions. The 
error estimates in |38| rely on classical results from Talay and Tubaro [36] and Milstein [28] 
based on the backward Kolmogorov equation (see also [291 Chap. 2.2, 2.3]). In this section, 
we denote by C^(M'^,M) the set of functions whose derivatives up to any order have 
a polynomial growth. We denote by C the infinitesimal generator of ( 1121 ) where we set 
fo = fi + f 2 - for any test function (j) G C^(M^,M), 


C(p 


a 


/o • + yAc)*, 


where V(/> denotes the gradient of (p and A(/> the Laplacian of (j). 


( 22 ) 


Theorem 2.6. fSSl Theorem4-1] Under Assumption\2^ let Xn he an ergodie numerical 
solution of m with hounded moments of any order M G N*, i.e. 

mXn\\^)< Cm (23) 


for all n > 0 , where Cm is independent of h,n. Assume further that the seheme has loeal 
weak order q > 1 i.e. it satisfies for all initial condition Aq = x and all h sufficiently small, 

\E{fiX,)) - E{f{X{hm < C{x, cP)h'i+\ (24) 

for all cj) G C^(]R'^,M), where x >->■ C{x,cj)) has a polynomial growth with respect to x. 

Assume also that the numerical solution admits a weak Taylor series expansion of the 
form 

E{4>{Xi)) = (j){x) + hAofix) + hfA\f{x) + ..., (25) 

for all (f G ,M), where At : C^(M^,M) —)• C^(M^,M), i = 0,1,2,... are linear dif¬ 

ferential operators with smooth coefficients. Let Gn denote independent and identically dis¬ 
tributed random maps inE^, independent of {Xj}j<n, with Xn = Gn{Xn) having hounded 
moments of any order, and satisfying a weak Taylor expansion of the form 

E{(l){Gn{x))) = (p{x) + h'^Aqfix) + C)(/i®+^), (26) 

for all cj) G C^(M'^,M), where the constant in O has a polynomial growth with respect to x. 
Assuming further 

iAq+[£,Aq]rp = 0, 

where we use the commutator notation [A, B] = AB — BA. Then the postprocessor Xn = 
Gn{Xn) yields an approximation of order r = q-\-l for the invariant measure. 


lim a. s. 

M^OO 


M 


1 ^ _ r 

+ 1 


(j){y)dp 

OO (y) 




(27) 
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< K{(l),x)e-^^^ +C{^)h'i+^ (28) 


E{(j){Xn)) - / (l){y)dfioo{y) 

for all 4> € with tn = nh, where X, K{(j),x),C{4>) are independent of n and h 

assumed small enough. 

Remark 2.7. We emphasize that the boundedness of moment condition (|23p can he easily 
proved for Runge-Kutta type methods, such as the proposed method d, following the 
methodology of (see also f2^ Chap. 2.2]) using the global Lipschitz continuity of the 
SDE fields. We also refer to where this assumption is discussed in the context of 
postprocessed integrators for SDEs. Notice also that the Lipschitz condition on the SDE 
fields and the ergodicity assumption on the numerical method could be relaxed using the 
concept of rejecting exploding trajectories, as introduced in m, see also | 2 f in the context 
of ergodic SDEs. 

We are now in position to state the main result of this section. We show that the new 
method (I13p satisfies the assumptions of Theorem 12.61 with q = 1, and thus has order r = 2 
of accuracy for the invariant measure of ergodic SDEs. 

Theorem 2.8. Under Assumption consider the method d with postprocessor Xn- 
Assume that Xn is ergodic when applied to the system d- Then Xn has order two of 
aceuraey for the invariant measure, preeisely, 

1 ^ _ r 

“ / (t>{yWoo{y) 

M^oo M + 1 ^ J^N 

E{4>{Xn)) - [ 4>{y)dpoo{y) 

JRN 

for all (j) G with tn = nh, where X, K{(j),x),C{(j)) are independent of n and h 

assumed small enough. 

The proof of Theorem 12.81 relies on the following lemma where conditions of order two 
of accuracy are derived for a perturbation for the linearized Euler method. The proof of 
this lemma is postponed to the Appendix. 

Lemma 2.9. Consider the following modification of the linearized Euler scheme for d, 

+ a2a\/hin) + {1 + a3hf[{Yn))aVhfn 

Yn = Yn+bihh{Yn)+b2hf2{Yn)+CaVh^n. (31) 

where ui, 02 , 03 , 6 i, 62 , c are fixed real eoeffieients. If the following conditions hold, 

T 2tti A bi — = 0, Ui A U 3 Y bi — = 0, a ‘2 A b2 — = 0, 

— - + 02 + 62 — = 0, (f >2 — ^l)[/ 2 ,/i] = 0, (32) 

then, assuming the ergodicity of Yn in d, the postprocessed scheme d satisfies the 
assumptions of Theorem \2 . 61 ivith q = I, and Yn has order two of accuracy for the invariant 
measure, i.e. it satisfies (l2^ . (l30P (with Xn replaced by Yn)- 


< C{f))h‘^ (29) 

< K{(P,x)e-^^- AC{(j))h\ (30) 
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Note that the Lie bracket [/ 2 , /i] = f^fi — / 1/2 involved in the second order conditions 
(I32p vanishes only when the flows associated to the fields /i, /2 commute, which is not true 
in general. We thus impose 61 = 62 - Still, the system (I32p has infinitely many solutions. 
Setting 61 = 62 = 0 for simplicity of the postprocessor, two solutions remain. Choosing the 
solution which minimizes the absolute value of ui and 03 , we obtain the following choice of 
coefficients for the order two scheme (j3ip . 


-2 + 

— —“3 — -- 


02 = C = 


1 

2 ’ 


61 = 62 = 0. 


(33) 


Proof of Theorem 12 . It is sufficient to prove that the method (I13h satisfies the same ex¬ 
pansions and (f 2 U]) with g = 1 as method (f5Tp . (l55p (with the same differential operators 

.Ao = T, Ai, and A.i). Indeed, applying Theorem 12.61 with g = 1 to the scheme (fT^ . we 
deduce that method (jl3l) also has second order of accuracy for the invariant measure, which 
concludes the proof of Theorem 12.81 

To recover the scheme (fT3P from ([3Tp , ([33]) , in the first line of (f3T]l one has to replace 
with Jn, 2 f.n = (d + 0{h))^ri hi the arguments of /i, /2 and also to substitute I + 03 / 1 /((x) 
with 


^l-V2 + V5 , 1 +V2-V5 

[ o -'n,l + r) 


)Jn,2 = I + azhf[{x) +0{h^); 


in the second line of m, one has to replace with Jn,‘ifn = {I + 0{h))f^n- We obtain 
that the difference between one step of (fT^ and one step of ([3Tp . (l33]) with initial condition 
Xq = yQ = X has the form Xi —Yi = R{x)f^hf’A _|_ 0{h^). Using E(^) = 0, we deduce 
E(<(.(Xi)) - E(<(.(yi)) = 0(/i3), while E(<(.(Xo)) - E((()(Fo)) = 0{h^). □ 


It can be seen from the proof of Theorem 12.81 that the operator J „^2 in front of and 
the operator J „^3 in the definition of the method (1131) have no influence on its order two of 
accuracy for the invariant measure in finite dimension. In infinite dimension, however, these 
operators play an important role for the well-posedness, the stability and the accuracy of 
the scheme in the SPDE case presented in Section [3l 

3 Analysis in the SPDE case 

3.1 Abstract setting and assumptions 

The state space in the SPDE case is an infinite dimensional separable Hilbert space TL, for 
which we denote by (•, •) the scalar product, and by | • | the associated norm. Consider the 
linear operator A involved in the parabolic SPDE ([1]). Recall that we assume that —A is 
an unbounded self-adjoint linear operator with eigenvalues 0 < Ai < ... < Ap < Ap+i < ..., 
such that Ap —)• -|-oo, when p —>■ -|-oo, and associated normalized eigenvectors ep (such that 
Aep = — ApCp), which form a complete orthonormal system in TL. 

For any s G M'*', we classically define the unbounded linear operator {—AyA from TL to 
TL and its domain TL^ G TL as follows: 

+CXD +CXD 

(—= ^(u,ep)Ap/^ep for all u E TL^ = {u G TL : |u|g = ^ |(u,ep)|^Ap <-boo}. 

p=l P=1 
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We also define the bounded linear operator {—A) and the semi-group both 

as linear operators from ^ to ^ by 

+00 +00 

{-Ay^'^u = '^{u, ep)X~'^/‘^ep, ^ exp(-tAp)(u, ep)ep. 

P=1 p=l 

It is straightforward that for any t £ (0, -|-oo) we have %) - the space of bounded 

linear operators from % to %, endowed with the norm denoted by || • || ~ with ||e*"^|| < 
exp(—Alt). Moreover, the following regularization property holds true; ||(— 
where Cg = sup^gR+ exp(— 2 r)r® £ ( 0 , -t- 00 ). 

Covariance operator We assume that Q is a bounded, non-negative self-adjoint linear 
operator from Ti to T-L, which satishes Qep = QpCp for any p G W, where the eigenvalues 
(^p)peN* ^ bounded sequence of non-negative real numbers. We assume that condition 
(j2]) is satisfied. The Q’Wiener process in ([1]) is then defined as follows: for any t > 0, 

+CXD 

W^{t) = J2V^Pp{t)ep, (34) 

p=i 

where (/3p)pgj^* is a sequence of independent standard scalar Wiener processes on an un¬ 
derlying probability space (n,T’, P). 

Notice that the operators A and Q commute: AQu = QAu for any u £ This 
property is often assumed in the literature, and simplifies the analysis of the order of 
convergence made in Section [T3l Nevertheless several arguments (especially ProDosition l3.ll 
and the results of Sectiondo not require this property; in particular, the scheme remains 
well-defined in the non-commuting case. 


Nonlinearity The nonlinear coefficient F is assumed to be a Lipschitz continuous func¬ 
tion from % to TL, with Lipschitz constant L satisfying the dissipation condition 


L = sup 


\F{u^) - F{u^)\ 
\u^ — v?\ 


< min Ab = Ai. 

pSN* 


(35) 


This condition ensures ergodicity of the continuous-time process lProDosition l3.ll) and of its 
time-discretized approximations (Proposition [TS]). A typical example of such a Lipschitz 
function on "H = L‘^{T>) - where F is an open smooth bounded domain in - is the 
Nemytskii operator F : u ^ f o u, where / : M —)• M is a globally Lipschitz function. 
Note that F = —DV is the derivative of the potential function V : L^{F) —)• M where 
V{u) = — fjy u{x) f {0u{x))dxd9. 

Under the above hypotheses, and assuming the trace condition ([2]), the process 
takes values in for any s < s, and we recall without proof the following result of expo¬ 
nential convergence to a unique invariant distribution, see e.g. [10] for general results, and 
[T3l Section 3.1.1]. 

Proposition 3.1. Assume (l2|) and the above hypotheses. Then, the process {u{t)'j 
solution of dH) admits a unique invariant probability distribution on TL. Moreover for 
all s < s, 

/ 1^^181^00(^1^^) < +00, 

JH 
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and for all 4> : Ti Lipschitz continuous, and all t > 0, 


E[(j){u{t))] - [ (j){v)pLoo{dv) 

Jn 


<C{(j),uo)e 


where C{4>,uo) is independent oft. 


Condition (I35h is crucial for the proof of the uniqueness of the numerical invariant dis¬ 
tributions established in the next section: we compare the solutions starting from different 
initial conditions and driven with the same noise process, and show an exponential con¬ 
traction similar to the result of Proposition 13.11 Notice that weaker conditions than (1351) 
are known in the literature (see e.g. m and references therein) to ensure the ergodicity 
and the exponential convergence of ([T|) - for instance when F is bounded and Lipschitz 
continuous with no size restriction on L. 


3.2 Stability and ergodicity of the integrator for SPDEs 

In this section, we prove the existence and uniqueness of invariant distributions for the 
time-discretized processes defined by the numerical method, for any time-step size h > 0, 
in a general setting. Notice that the results of this section do not require the gradient 
assumption F = —DV. The results are analogous to classical results for the 0-method in 
the context of stiff SDEs [HI Theorem 3.1] and for the linearized implicit Euler method ([^ 
in the context of SPDEs as studied e.g. in [U Remark 4.8]. 

It is a key observation here to exploit that the sequence (u„,Mn-i)neN defining the new 
scheme dSD,® is a Markov chain on the product space T-L x T-L. The initial condition 
{uo,U-i) is given by uq = tt( 0 ) = uq and an arbitrary It_i G H which plays no role in the 
dynamics, since {un+i,Un) depends only on Un and ^n, not on Un-i- 

In the following proposition, we state uniform bounds - with respect to n G N and 
h G (0,1) - on first-order moments for the norm j • jg for s < s. 

Proposition 3.2. Assume the hypotheses of Section \3.1\ and consider the scheme ©,(110]). 
For all s G [0,s), assuming uq G there exists a constant Cs G (0,-|-oo) such that for all 
h G (0,1), 

supE|u„|^ < (73(1 -F Juols), supE|m„|^ < (73(1 -b luojs). 

nSN nSN 

Proof. Thanks to ([21), Trace(JjQJi) < -boo for i G {1,2,3}, and thus Un and Un are well- 
defined in PL for all n G N. The contributions of the drift part and of the stochastic 
perturbation are treated separately: we introduce the auxiliary process as the 

solution of the following equation 

4+1 = + cTVh{^^^^Ji + ^ 7 ^/) 

with 4 = 0. Set dn = Un — 4 for all n G N; then do = uq and 

4+1 — Jldn “b hJ^Fi^dn “b 4 “b 
The quantity in satisfies for all n G N* the identity 

in — crvhy 2 7i -b ^ 4' 

k=0 
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Observe that ^ ^ is a bounded linear operator from "H® to 'W with 

norm less than ^ for s G [0,s). Since the Gaussian random variables {^k)k£N 
independent, for all 0 < s < s, 


^ ^ k=0 


n—k^Q\2 

fc Is — 




QpK 


^ 2 +00 


< 


Xph 


3- 


3 -v 2 ,Gi,tr(i+w‘ 

Trace((—. 


< 


(1 + Xph)"^ — 1 3 — y/2 


P=i ' -'P' 

Now thanks to (j35l) . straightforward computations show that 

]E|d„| < + :rTWl^(0)l + T^^{Wn-l\+E\lj2CTVhCt,\) 


< 


1 + Xih 
(1 + hLf' 
(1 + X\hy 


1 + X\h 

>o| + + 


L 


Ai - L Ai - L 


1 + X\h 
(supIE|4| + 


1 


fj(Trace ((-A) 


km 


3 - ^/2 


with < 1 — h < exD < 1 

(l+Ai/i) - l+Aifc ^ -‘-• 

As a consequence the claim follows for {un)nm iii the case s = 0. In particular, for some 

constant C G (0,+oo) it comes that sup^gpjE|F(un + ^J 2 (^Vh^n)\ < C{1 + |uo|). 

The case s G (0,s) is treated using the estimates of Lemma 3.2 in [3]: for n G N 

n—1 - 

k=0 

- (1 +A,+ '>^'(1 + “»!) S( (^p5 + ’‘“>‘(i + a,A)‘-W/.j) 

< (78(1 + IuqIs)- 

Finally, for n G N, 

= ^Trace((/- ^A)-i(-A)^Q) < ^Trace((-A)-i+^Q). 

This concludes the proof, since Un = Un + ^J^crVh^n ■ D 


We now state the following existence and uniqueness result for the invariant distribution 
of the Markov chain {^Un,Un-i)j^^^- We prove that the convergence is exponentially fast, in 
contrast to the trapezoidal method (|lip which is not L-stable and for which such exponential 
estimate does not hold in general for stiff problems (even in finite dimension). 

Proposition 3.3. Assume the hypotheses of Section [Ql For any h G (0,1), the Fi x FL- 
valued Markov chain (?x,i, admits a unique invariant distribution inFL xFL, 

with marginals in FL denoted by and respectively. 

Moreover, the convergence of the distributions to equilibrium is exponentially fast: for 
all Lipschitz test function ip : FL m M., and all tn = nh, 


Eip{un) - / + E^iun) - / FdFx’fo < l^^ol) exp 

Jh Jh 


(Ai-T) ^ \ 

1 + Ai/i "y ■ 
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Proof. Existence. The semi-group on PL x PL generated by the Markov chain 

(ri„, satisfies the Feller property: for any n £ N, for any bounded continuous test 

function (/) : x —>• M, the map (uoi^-i) P^^f>{uQ-,U-i) = lE(/>(un, u^-i) is continuous. 
The claim then follows from the standard Krylov-Bogoliubov criterion (see Section 3.1 in 
[lOjl: given an arbitrary initial condition {uq,U-i) £ PL x PL, if M„ denotes the law of 
{un,Un-i), then 

• (^ELoMi)„a IS a tiQht seQuence of probability distributions on PL x PL as 
a consequence of Proposition 13.21 combined with the Markov inequality, and of the 
observation that for any s £ (0,s) and any i? > 0, the set {|u|s < R, |uls < R} is a 
compact subset of PL x PL. 

• every subsequence limit point M is an invariant distribution for the semi-group. 

Uniqueness. Consider two initial conditions UqjUq £ PL, as well as £ PL and 

the associated processes and * = 1)2, defined by ®, cni), and driven 

by a unique noise process (Cn 

Then by Lipschitz continuity of F, and using the cancellations of several noise terms, 
computations similar to those of the proof of Proposition 13. 2l vield for any n £ N the almost 
sure contraction property 


Wi - ul\ 


= \ui - ul\ 


< 


1 + Lh 
1 -|- Xih 


<exp - 


(Ai - L) 

1 -|- Xih 


Un - Un 


Finally, taking {uQ,ufi) random, independent of the noise process and distributed 

according to an ergodic invariant distribution M^ gives the exponential convergence and 
the uniqueness properties. □ 


3.3 Analysis of the order of convergence: a simplified linear case 

It is shown in [3] (for d = 1 and s = 1/2, associated with space-time white noise Q = I), that 
the standard linearized implicit Euler scheme @ has order r = 1/2 — e for all e £ (0,1/2) 
for the approximation of the invariant distribution /ioo of ([T]). In this section, we show 
that the postprocessed scheme has the improved order of convergence s -|- 1 — e. Since the 
techniques from Section [2] do not extend straightforwardly to the SPDE case, we only focus 
on a simplified case, where the nonlinear coefficient F is replaced with a bounded linear 
operator. Numerical experiments of Section |3] show that the higher order is preserved for 
various examples of nonlinearities F. 

In addition to the hypotheses of Section 13.11 assume that the coefficient F is given by 
a linear mapping: for any u £PL, F{u) = Bu where B £ C{PL) satisfies Bcp = —bpCp for all 
p £ N*, with real eigenvalues bp £ (—Ai, Ai) (due to condition ([551) 1: 

du{t) = Au{t)dt + Bu{t)dt + adW^ {t) , n(0) = uq. (36) 

In this situation, the components {u{t),ep) in the basis {cpjpgN* of the solution 
of the SPDE ([T|) are independent Ornstein-Uhlenbeck processes. Similarly, the components 
in the basis {epjpgN* of the discrete-time processes (resp. resp. 

are also independent processes. 
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As a consequence, explicit expressions for the invariant distributions and 

are available: they are centered Gaussian probability measures on Ti, 

2 

IJ-OO = AA(0,(3oo) , Qoo = yQ(-^ - 

= Af(0,Qy , Tt^-^{0,Ql). 

The expressions for and being more complicated are displayed in (ITO below. 

The following lemma is a key elementary tool in order to exhibit the order of convergence 
as /i —)• 0 of the approximating measures towards /Too- 

Lemma 3.4. Let TTj = M{0,Qj), j € {1,2} be two Gaussian probability distributions on 
the Hilbert space H. Assume that for all p £W, j G {1, 2} QjCp = qj,pep. Let ip G C^iH, M) 
satisfy that sup„g-^ ||Zl^(/9(u)|| < +oo. Then 


/ ipd7r2 — / ipd-Ki 

Jn Jn 


< 


sup^ew \\D‘^t{v)\\ 


+ 00 


1 ^ 2 ,P 

p=l 


9l,p| 


(37) 


Remark 3.5. Assume that 52 ,p > qi,p- Then the above Lemma \S.4\ yields optimal orders 
of convergence. Indeed choosing the test function ‘Popt{u) = exp(—|up), Lemma 9.5 in f2^ 
yields that 


Jh Jh 


Trace (<52 - Qi) 
exp(6Trace((52)) 


This means that the quantity Trace(Q 2 ~ Qi) o,lso provides a lower bound for the error 
between the invariant distributions vri and 712 - 


Proof of Lemma \3.4[ Let (7p)pgpj* and ((5p) be two independent sequences of i.i.d. 
standard real valued Gaussian random varia&es, centered and with variance 1. Set Xj = 

EpeN* = EpeN* ^max((-l)l(gi,p - q2,p),0)Spep, for j G (1, 2}. Observe 

that Xj ~ TTj, and that Xi + Ri and X 2 + R 2 have the same Gaussian distribution. This 
yields 


f (fdTT2- [ (fdTTl = E[(^(X2)] - E[(/?(Xi)] 

Jn Jn 

K[ip{X2)] -E[ip{X2 + R2)] +E[ip{X^ + R,)] -E[ip{X^)] 
Ey{X2 + R2)] -Ey{X2)] -E[Dip{X2).R2] 


< 


+ E 


[(^(Xi+Ri)] -E[yp(A:i)] -E[DgT{Xi).Ri] 


Indeed, E[Dip{Xj).Rj~\ = 0, because Xj and Rj are independent and E[i?j] = 0. Using a 
second-order Taylor expansion, we deduce 


Jn Jn ^ 


where we note E|i?i|^ -|- E|i? 2 |^ = Ylp^i\Q 2 ,p - qi,p\- 


□ 
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We now explain how Lemma 13.41 permits to find the order of convergence of to 
fioo- We define, for p G N*, the component processes u.{p) and u.{p) by projecting on the 
eigenvector e^: for any n £ N, Un{p) = {un, Bp), Un-i{p) = Bp). Then (P |) . (fTUl) applied 

to (IHHIi rewrites as a system of independent equations, decoupled with respect to p G N*, 


^n+iijp') — “ 4 ( bph^Uyii^p') + (j'\/hyJ Qpi3(^ 

Unip) — C{ ^p^)'Wn(p) T ^ph)^n,pj 


(38) 


nSN.peN* 


are independent standard Gaussian ran- 


where ^^n,p = {^n,Bp): thus {^n,p) 
dom variables, and the rational functions A,I3,C,'D satisfy for any z G (—oo,0) and 
/3 G (-l,min(l,| 2 |)). 


A{z,P) = ^ B{z,l3) = 


1 + - — 


/3 3-V2, 


(1-Z)(1- 


3-V2. 


Ciz) = 1, V{z) = 


2{1-zliyA' 

(39) 

Since for /3 G (—1, min(l, | 2 ;|)) the stability condition \A{z,j3)\ < 1 is satisfied, straightfor¬ 
ward computations yield 


Q^Bp= lim W\Un{p)\'^Bp = ^ lZ{-\ph,-bph) 


n—^+c» 


2 (Ap -|- bp) 


Q'Lep= lim E\un{p)\^Bp = 

n^+oo 2(Ap -|- bp 


— — — Tl{—Xph, —bph)Bp 


-2{z+p)B(z,P)‘^ 


where 7^(z, /3) = 


W 


and 


n{z, /3) = C{z)‘^n{z, /3) - 2{z + (3)V{z)‘^ = 1 + I3z 


Pyz)/3 + P2{z) 
{2 + /3-z)Pyz) 


(40) 


with polynomial functions Pi{z) = 10 — 4\/2 — (11 — 6y/2)z, ^ 2 ( 2 ;) = 20 — 8\/2 — (44 — 
24 ^ 2 ) 2 ; + (11 - 6V2)z^ and_P3{z) = (2 - z)(2 - (3 - V2)z)‘^. 

The following estimate P{z, j3) = l + 0{zj3) as z, ^ 0 is crucial to obtain an improved 
order for the convergence of to poo when /i —)■ 0. It is not surprising because the scheme 
samples exactly the invariant measure of (|36|) in both cases ^ = 0 or i? = 0 (as already 
shown in Proposition 12.311 . equivalently 7 ^( 2 ;, 0) = 7^(0,/3) = 1 for all z,/^. 

Lemma 3.6. For all z <0 and all (5 G (—l,min(l, |^|)), wb havB 


1 - 7^(2:,/3)| < \zP\ 


(15 - 6^2) 
4(1-2)2 • 


Proof. Observe that for all 2 < 0, we have P 2 {z) > Pi{z) > 0 and (2 + /3 — 2 ) > 0, 7 ^ 3 ( 2 ) > 0. 
Since (2 -|- ,5 — 2 )“^ < (1 — - 2 )“^, we obtain 


1 - 7^(2,/3)| < |2/3| 


Pi( 2 ) + P 2 (^) 
{l-z)P3{z) ■ 


The estimate then follows by observing that (1 — is an increasing function of 

2 < 0, with maximum at 2 = 0, given by = (i 5 - 6 \/ 2 ) ^ 1 ^ 
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We are in position to state our main convergence result, which yields the order of 
convergence r = s + 1 for any s < s for the invariant distribution with postprocessing . 

Theorem 3.7. Consider the method (l^- (fT0]l applied to (l36|) . Let ip £ C‘^{T-L,'R), such 
that sup^g-^ ||Z)^(/9(u)|| < + 00 . For any s < s there exists Cs £ (0,+oo) such that for any 
h G {0,1/L) we have 




< Cs sup \\D‘^(p{v)\\a'^h^^^. 
v(^n 


(41) 


Moreover, if B = 0 (i.e. bp = 0 for all p £ N*j then the method is exact: 
Proof. Thanks to Lemma 13.41 it is sufficient to control 


+ CXD +00 2 

|((Qoo - Q^)ep,ep)\ = ^ ^ 

p=i p=i 


+00 2 

\n{-Xph,-bph)-l\ 


H-cx) 


< CAicj 2 ^ 


< CXla^ ^ qpXr^+^ 


^ Xp {1 + Xphy 
< CsXla^TXace{{-A)-^+^Q)h^+^, 


-i+s,;_ 

^ (l + Xph)^ 


which gives the order of convergence s +1 for all s < s. Moreover, it is clear that TZ{z, 0) = 1 
for any z < 0, so that if 6p = 0 for all p £ N* then = Qoo- D 


Remark 3.8. We see in Lemma \3.6\ that the error is zero if z = 0 or (5 = 0. This is related 
to Proposition IS.,?! which shows that the error of the postprocessed method is zero in the 
linear case when fd = 0. This feature permits to gain one power of h and thus one order of 
accuracy in the proof above. In contrast, notice that the standard linearized implicit Euler 
scheme has the lower order s for all s < s. Indeed, under the hypotheses of Theorem [13 
then for all h small enough, Qoo ~ QSo u nonnegative and using Lemma \3lf\ we only need 
to control 


2 + 0 O 


■nace(Q„-*.J=TE 


2 ^^ Xp -j- 2 -j- Xph -\- bph 






-l+S 


A. 


’-h^ 


p=i 


Xn H” 2 -j- Xxyh -\- bnh 


< ^^C'sTrace((-A) 


Thus for all s < s that there exists Cg £ (0, +oo) such that for all h small enough 

/ ipdpoo - / < Cg sup ||L>^(/7(u)||(T^/iL 

Jn Ju ven 


(42) 


It is also possible to prove that the above Theorem \3. 7| (resp. (|4^ ) gives optimal order of 
convergence, namely that m does not hold true for all h > 0 if s > s + 1 (resp. (I42|l 
does not hold true for any h > 0 if s > s). This fact is also supported by the numerical 
simulations of Figure 0 in Section 
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3.4 Spatial regularity analysis 


We show in this section that the postprocessed method yields a solution which has the same 
regularity in space as the exact solution, in contrast to the standard linearized implicit Euler 
method, which yields a solution that is too smooth. The action of the postprocessing thus 
not only increases the order of the convergence, but also provides a qualitatively better 
approximation with the correct regularity. For all Borel probability measure fj, on T-L, we 
define its regularity, denoted reg(^) € MU {—oo. Too}, by the supremum of s such that the 
norm | • |s of is square-integrable with respect to /r: 

reg(/r) = sup{s G M, / < oo}. (43) 

Jn 

The interpretation in terms of random variables is the following. For a random vari¬ 
able V with values in "H, denoting its probability law, we have the identity E(|ul^) = 
|u|gP.u((itt) which is a finite quantity if s < reg(P.u) and -|-oo if s > reg(P.i,). Notice that 
instead of quantifying the regularity in terms of the Sobolev space one could state sim¬ 
ilar results in terms of Holder regularity. We refer to [9] for such a study of the 0-method 
applied to the stochastic heat equation with finite differences. 

We focus for simplicity on the case F = 0, a = 1, and with the initial condition uq = 0, 
but we emphasize that the extension to the general semilinear situation is straightforward, 
the regularity being determined only by the stochastic terms in our setting. First, notice 
that for the exact solution u{t), the regularity parameter is reg(Pu(t)) = reg(/roo) = s for 
all t > 0, (see also Proposition 13.ip . Indeed, using ([3]), and the Ito formula yields 


E|«Wls 


+00 

E 


p=i 


Qp 


(1 


— exp(—2Apt)) 


< -1-00 
= -|-oo 


if s < s 
if s > s 


The following proposition shows that at the discrete-time level, the standard linearized 
implicit Euler method Vn in Q and the method without processing Un in Q have the 
regularity parameter s-|- 1, whereas the postprocessor Un in (Hop has the correct regularity 
parameter s. 


Proposition 3.9. Consider ([ip with F = 0, cr = 1, uq = 0 and assume djp. Then, for 
all h > 0 and all n G N*, reg(P„„) = reg(Pu„) = reg(i/^) = reg(/r^) = s -|- 1, whereas 
reg(Pji„) = reg(F^) = s. 


Proof. Inspecting the proof of Proposition 13.21 we have Un = and for all s < s -|- 1, h > 
0,n G N* 




- 1-00 


(Ap/i)" 




+ - lTrace((-.4)-Q); 


this yields reg(Pu„) = reg(/i^) > s -|- 1. The reverse inequality is obtained with a similar 
lower bound. The proof for the standard linearized implicit Euler scheme Vn is similar. 
Now, adding the postprocessing, Un = Un + ^J^cry/hf^n , we obtain using the definition of 
^3, 




- 1-00 


hTtace{{I - -A)-\-ArQ) =Y, 

p=i 


Qp 


Xph 


< - 1-00 


Ap ^ 1 -j- ^Xph ] — “t-oo 


if s < s, 
if s > s. 


The term ^J^aVhfn thus has exactly the same regularity s as the exact solution. This 
concludes the proof of reg(Pp^) = Teg{]l!^) = s. □ 
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4 Numerical experiments 


In this section, we compare numerically the performances of the new postprocessed method 
(1131) with the standard linearized implicit Euler method ©, and the trapezoidal method 
m, both in finite and infinite dimensions. 

We consider first in Figure [T] the scalar nonlinear SDE (|14|) with dimension = 1, 
parameters A = —l,cr = 1 and the initial condition X(0) = 0. Taking fi{x) = Ax and 
f 2 {x) = f{x) in (fT3l) . we consider the nonlinearities f{x) = —x—sin(x) and f{x) = —2x—x^, 
respectively, and we compute the averages over 10 ^® independent trajectories with final 
time T = 1 and compare for many time stepsizes the accuracy for E(exp(—X(r)^)) = 
exp(—The final time T = 1 is chosen large enough so that the equilibrium 
is reached and the exponentially decaying term in (|3np is negligible. In the left 

picture of Figure [T] where the nonlinearity f{x) = —x — sin(a;) is Lipschitz, we observe 
as shown in Theorem 12.81 the expected order 2 of convergence for the new method, while 
the standard methods exhibit order 1 of convergence (see the reference lines with slopes 
1,2). Although our analysis in Section f2.2l applies only to globally-Lipschitz vector fields, 
we observe that the excellent performances of the new method persist also in the example 
with the non-Lipschitz nonlinearity f{x) = —2x — x‘^ and the globally bounded test function 
(/)(x) = exp(—(right picture of Figured]). 

We next consider a standard finite-difference approximation Uj{t) ~ u{jAx,t) of the 
ID heat equation ([d|) with zero Dirichlet boundary conditions on a uniform grid with size 
Ax = 1/(A^ -l- 1). This yields the following system of SDEs in dimension N, 



(-2 1 \ 






fdW\t)\ 

1 

1 -2 1 


X^{t) 


fixHt)) 

, 1 

dW\t) 





dt “h 


dt -|- - 

V Ax 



V 1 -V 


[x^{t)) 




\dW^{t)) 


where W ^,..., are independent one-dimensional standard Wiener processes. 

Considering N = 100 grid points for the space discretization, we take the initial condi¬ 
tion rt(x,0) = sin(27rx) and plot in Figure [2] a sample trajectory on the time interval (0,1) 


Case f{x) = —X — sinfx) Case f{x) = —2x — x^ 



Figure 1: Comparison of the new method (solid lines) with the standard linearized 

implicit Euler method (dashed lines) and the trapezoidal method (dashed-dotted lines) 
for the scalar SDE p4|) (dimension N = 1) with A = —1,(7 = 1, nonlinearity /(x). 
Error for E(exp(—A(T)^)) at final time T = 1 versus the stepsize h, where 1/h = 
8,12,16, 24, 32, 44, 64,92,128. Averages over 10^° samples. 
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(a) Standard linearized Euler method, sample trajectory u(x,t) and corresponding profile at 
final time t = 1. 



(b) New method with postprocessor, sample trajectory u{x, t) and corresponding profile at final 
time t = 1. 


Figure 2 : Samples of realisation of the stochastic nonlinear heat equation dH) with f{u) = 
—u — sin(u) using the standard linearized implicit Euler method and the new method with 
postprocessor. With N = 100 space grid points and timestep size h = 1/100. 


for the standard linearized implicit Euler method ([ 6 |) and the new method with postproces¬ 
sor Xn, using the same sets of generated random numbers for both methods. We observe 
that the solution of the standard linearized implicit Euler method (Fig.[2]^a)) is qualitatively 
too smooth compared to the new method (Fig.[2]^b)) with the postprocessor applied at each 
timestep, which corroborates the statement of Proposition 13.91 in Section 13.41 Notice that 
the trajectory for the new method without applying the postprocessor would look very 
similar to that of the standard linearized implicit Euler method. The spatial regularity 
observed in Eigure [2]^b) is the same as the one of a diffusion process driven by Brownian 
motion, conditioned to be zero at initial and final times. This property is not surprising, 
see dalle] for a study of the link between the law of a conditioned diffusion and invariant 
distributions of SPDEs. 

We finally plot in Figure[3]the error at final time T = 1 for the quantity E(exp(—||u||^ 2 (o i))) 
where we use the approximation 11 ^* 11 ^ 2(0 i) — arbitrarily the initial 

condition u{x, 0) = 0 and we compute the averages over 10® independent trajectories so 
that the Monte-Carlo errors become negligible. The reference solution is computed using 
the new method with stepsize h = 1/512. We consider respectively, the cases of various 
nonlinearities f{u) = 0, /(«) = —u, /(«) = —u — sin(u), f{u) = —2u — v?. We observe 
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Case f{u) = 0 Case f{u) = —u 



Figure 3; Comparison of the new method (solid lines) with the standard linearized implicit 
Euler method (dashed lines) and the trapezoidal method (dashed-dotted lines) for the 
stochastic heat equation (j3|) with nonlinearity f{u) discretized in space with N = 100 
grid points. Error for E(exp(—||rt(T)|p)) at final time T = 1 versus the stepsize h, where 
ljh = 8,12,16,24,32,44,64,92,128. Averages over 10® samples. 


in all cases an order of convergence 1/2 for the standard linearized implicit Euler method, 
while the trapezoidal method has order 1. In contrast to the standard linearized implicit 
Euler method, the new method with postprocessor has a much better accuracy by a factor 
15 — 250 in the range of stepsizes considered. It has a zero bias for f{u) = 0 and order 
3/2 in the linear case f{u) = —u, as proved in Theorem 13.71 We observe that the order of 
convergence persists in the nonlinear case f{u) = —u — sin(ri), and the excellent accuracy 
persists in the non-Lipschitz case f{u) = —2u — u^, although an order reduction can be 
observed. 

Acknowledgements. The research of C.-E. B. and G. V. is partially supported by the 
Swiss National Science Foundation, Grant: No 200020-149871/1 and No 200020_144313/1, 
respectively. The computations were performed at University of Geneva on the Baobab 
cluster. 


Appendix 

Proof of Lemma \2.!A We consider the variant (1311) of the linearized implicit Euler method. 
A straightforward calculation yields the following weak expansion for all (f £ (M^, M), 

E((;i)(Yi)|lo = x) = 4>{x) + hCcj){x) + hfAict){x) + 0{h?), (44) 
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where the constant symbolized by O is independent of h but depends on (j) and depends on 
X with a polynomial growth. Here, the fourth-order linear differential operator Ai is given 

bjQ 

2 N 4 Af 

2=1 'i'd=^ 

2 N N 

+ ^(ai + 1)^^ —4>' ^ Gi) + + 1 + 03 ) 0 -^ ^ Cj) 


2 = 1 


2=1 


N 


2 N 

2=1 2=1 


where ei,..., ew denotes the canonical basis of R-^. 

For any V’ G C^(R^,R), we denote (ip) = J^n P^{y)f^oo{dy) = J^n ^P{y)p{y)dy. By 
integration by parts, and using the identity Vp = -^pfo, we obtain, 

(</'"(/o,/o)) = (- 0 '(/o/o + (div/o)/o + ^ll/olP/o)), 

Ei ei, Ci)) = {p'ia^ Yli foi^i, Ci) + 4/^/o + 2 (div /o)/o + ^||/o|P/o)) , 

= (- 2 (/>'"(/o, e^, e*)) , 

</>"(/(e^ei)) = (-p'(<^^I^ifi(ei,ei) + 2flfo)), (45) 

see mm for examples of such calculations. Then 

/ 2 ^ 1 ^ 

= ((a?+ 2ai^y(()'^/"(ei,ei)-b (^-+ ai+a3^o-2^0"(/(ei,ei) 

\ i=l i=l 


+ alY(p''^fUei,ei) + ( - 1 + 02)0-^ ^ Ci) 


N 


N 


2=1 


2=1 


We note that the postprocessor Yn Yn in (fdTTi satisfies (l26li with q = I and Aip = 
hP'fi + b2(t>' f2 + wO-2A(/). 


N 


N 


[C,Ai\(t) = {bi - b2)(j)'{f2fi - f[f2) - —o-2(/>'^/Q'(ei,ei)-c2o-2^(/)"(/^ei,ei). 

2=1 2=1 

We deduce the following expression for (^Aip + [£, ^i]^), which is the key quantity to 
compute in Theorem 12.61 Summing up, we obtain 

{Alp A [C,Ai\p) 

2 N N 

= (^(^al + 2ai +bi- ^/"(g-, gj) + (ai + as+ - +bi - ^ (/>"(/(e*, e*) 


2=1 


2=1 


2 N N 

+ ^02 + 62 ~ ~ 2 ‘^' E (^*5 ^i) A (^ — — A a 2 A b 2 — ep) 

i=l i=l 


^ We denote 4>'{x) : R the first derivative of <j} af point x € R^, 4>"{x) the second derivative (a 

symmetric bilinear form on R-^ x R^), 4i'"{x) the third derivative (a symmetric bilinear form), etc. 
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The above quantity is zero by assumption for all test function cj). This means (^i + 
\C,Ai])*p = 0. Applying Theorem 12.61 with g = 1 to the scheme (I^TI) then yields for Yn 
in (IdID a method of weak order two for the invariant measure. This concludes the proof of 
Lemma 12.91 □ 
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